Acquisition

Author

Claude Grasland

Le but de ce chapitre est d’approfondir les méthodes de recueil de données lorsque celles-ci comportent une information géographique sur la localisation des objets étudiés. Nous allons donc charger non plus seulement des tableaux statistiques mais aussi de l’information géographique décrivant la localisation de poins, lignes ou polygones. Cela implique deux nouveautés :

Nous verrons en détail l’utilisation des données géographiques dans les chapitres ultérieurs et on se bornera ici à analyser comment recueillir cette information.

PREPARATION DU TRAVAIL

Packages

On charge les packages habituels

knitr::opts_chunk$set(echo = TRUE, warning= FALSE, message = FALSE, error=FALSE)

## Affichage de tableaux
library(knitr)

## Requêtes web
library(httr)
library(jsonlite)

## Tidyverse & co
library(dplyr, warn.conflicts = T, quietly = T)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(ggplot2)

## Information géographique
library(geojsonsf)
library(sf)
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

Organisation du projet

On crée un dossier resul à l’intérieur duquel on va imporer les données statistiques (sous-dossier stats) et des données géométriques (sous-dossier geom). On ajoute un dossier sf où on placera les ficheirs au format spatial features de R.

Objectif Grand Paris

On se donne comme ojectif d’acquérir des données géométriques décrivant les communes et les iris du Grand Paris qui regroupe quatre départements 75, 92, 93, 94 et quelques communes isolées de la grande couronne dont les codes sont : “95018”,“91027”,“91479”,“91432”,“91589”,“91326”,“91687”.

Nous devons arriver à produire quatre fonds de carte :

  • par IRIS
  • par communes
  • par territoires
  • par département

Il faudra faire en sorte que les résultats se superposent …

Importation de fichiers SIG

Beaucoup de données géographiques ont été conçues pour l’utilisation de systèmes d’information géographiques tels que ArcGIS ou QGIS. Historiquement ces données sont le plus souvent enregistrées et échangées dans le format shapefile qui se compose de trois ou quatre fichiers pour chaque fonds de carte :

  • carte.shp : fichier contenant la géométrie (contour des unités)
  • carte.shx : index des unités géométriques
  • carte.prj : fichier contenant la projection de la carte
  • carte.dbf : fichier de données attributaires (statistiques)

Nous allons prendre à titre d’exemple l’acquisition d’un fichier des communes d’Ile de France disponible sur le site data.gouv.fr à cette adresse

Nous commeçons par télécharger le fichier avec la fonction download.file() et on stocke le résultat dans notre dossier “geom”.

myurl <- "https://www.data.gouv.fr/fr/datasets/r/5cd27d86-4859-40dc-b029-a215219eedf9"
download.file(url = myurl, destfile = "resul/geom/idf_com.zip")
list.files("resul/geom")
[1] "GP_com.geojson" "GP_com.RDS"     "idf_com.zip"   

Nous avons téléchargé un fichier au format .zip que l’on va décompresse à l’aide de la fonction unzip()

zipF<- "resul/geom/idf_com.zip"
outDir<-"resul/geom/"
unzip(zipF,exdir=outDir)
list.files("resul/geom")
[1] "communes-dile-de-france-au-01-janvier.dbf"
[2] "communes-dile-de-france-au-01-janvier.prj"
[3] "communes-dile-de-france-au-01-janvier.shp"
[4] "communes-dile-de-france-au-01-janvier.shx"
[5] "GP_com.geojson"                           
[6] "GP_com.RDS"                               
[7] "idf_com.zip"                              

On voit bien apparaître les quatre fichiers qui définissent un shapefile et on pourrait les utiliser avec des logiciels tels que ArcGis, Qgis ou Magrit.

Importantion de shapefile -> R

Nous allons maintenant importer le fonds de carte dans R à l’aide du package sf (spatial features) en utilisant la fonction st_read() :

map<-st_read("resul/geom/communes-dile-de-france-au-01-janvier.shp")
Reading layer `communes-dile-de-france-au-01-janvier' from data source 
  `/Users/claudegrasland1/git/datamining2024/resul/geom/communes-dile-de-france-au-01-janvier.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1287 features and 9 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 1.44728 ymin: 48.12054 xmax: 3.555687 ymax: 49.23719
Geodetic CRS:  WGS 84

On vérifie la classe de l’objet :

class(map)
[1] "sf"         "data.frame"

On voit qu’il s’agit à la fois d’un data.frame et d’un objet de type sf. Regardons plus en détail avec str()

str(map)
Classes 'sf' and 'data.frame':  1287 obs. of  10 variables:
 $ objectid  : num  12 14 56 57 59 73 78 86 90 94 ...
 $ shape_leng: num  16585 17790 14908 6894 10681 ...
 $ insee     : num  77472 91315 78472 77038 78672 ...
 $ nomcom    : chr  "La Trétoire" "Itteville" "Orsonville" "Boissettes" ...
 $ numdep    : num  77 91 78 77 78 95 95 95 92 77 ...
 $ fusioinsee: chr  NA NA NA NA ...
 $ nomcomto  : chr  "Trétoire (la)" "Itteville" "Orsonville" "Boissettes" ...
 $ st_areasha: num  9488102 12122472 9504776 1603038 5117345 ...
 $ st_lengths: num  16585 17790 14908 6894 10681 ...
 $ geometry  :sfc_POLYGON of length 1287; first list element: List of 1
  ..$ : num [1:17, 1:2] 3.26 3.26 3.24 3.24 3.23 ...
  ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA
  ..- attr(*, "names")= chr [1:9] "objectid" "shape_leng" "insee" "nomcom" ...

Par rapport à un dataframe classique on trouve une nouvelle colonne appelée “geometry” ainsi que différents attributs tels que la projection. Nous trouvons donc en un seul fichier l’ensemble des informations normalement présents dans les quatre fichiers qui composaient le shappefile.

Visualisation de la geometrie

Pour visualiser rapidement le fonds de carte, il suffit de taper la fonction rbase plot appliquée à la colonne geometry du fichier sf :

plot(map$geometry)

Extraction du data.frame sans les données géométriques

Si l’on veut juste travailler sur les données statistiques du tableau, on peut retirer la géométrie du fichier pour revenir à un pur objet de type data.frame à l’aide de la fonction sf st_drop_geometry() :

don<-st_drop_geometry(map)
class(don)
[1] "data.frame"
head(don)
  objectid shape_leng insee               nomcom numdep fusioinsee
1       12  16584.790 77472          La Trétoire     77       <NA>
2       14  17789.557 91315            Itteville     91       <NA>
3       56  14908.000 78472           Orsonville     78       <NA>
4       57   6894.431 77038           Boissettes     77       <NA>
5       59  10681.145 78672  Villennes-sur-Seine     78       <NA>
6       73  18855.711 95541 Saint-Clair-sur-Epte     95       <NA>
              nomcomto st_areasha st_lengths
1        Trétoire (la)    9488102  16584.790
2            Itteville   12122472  17789.557
3           Orsonville    9504776  14908.000
4           Boissettes    1603038   6894.431
5  Villennes-sur-Seine    5117345  10681.145
6 Saint-Clair-sur-Epte   12431151  18855.711

On constate que parmi les colonnes il existe à la fois une variable relative au département et une autre au code INSEE des communes. On peut donc procéder à une sélection sur ces deux critères afin d’aboutir à un fonds de carte du grand Paris.

Extraction des communes du grand Paris

On utilise la condition “ou” pour filter à la fois sur les départements et les communes.

GP_com <- map %>% filter(numdep %in% c(75, 92, 93, 94) |  
                           insee %in% c(95018,91027,91479,91432,91589,91326,91687))
plot(GP_com$geometry)

Cartographie rapide du résultat

On utilise une petite astuce pour visualiser les départements. On commence par créer une variable de type factor et on applique un plot spécial sur cette variable.

GP_com$Dep<-as.factor(GP_com$numdep)
plot(GP_com["Dep"])

Sauvegarde du fichier sf

On enregistre le fichier dans le format interne de R afin de ne pas avoir à répliquer toutes les étapes précédentes à l’aide de la fonction saveRDS(). On en profite pour nettoyer le dossier en supprimant les fichiers dont on n'a plus besoin à l'aide de la fonctionfile.remove()pour un fichier unique ouunlink()`` pour un groupe de fichiers

file.remove("resul/geom/idf_com.zip")
[1] TRUE
unlink("resul/geom/communes*")
saveRDS(GP_com,"resul/geom/GP_com.RDS")

Exportation au format Geojson

Si l’on doit travailler avec des programmeurs qui utilisent Python, le plus simple est de réaliser une exportation au format geojson à l’aide de la fonction st_write() du package sf. On rajoute l’option delete_dsn=T pour écraser une éventuelle version antérieure.

st_write(GP_com, "resul/geom/GP_com.geojson",delete_dsn = T)
Deleting source `resul/geom/GP_com.geojson' using driver `GeoJSON'
Writing layer `GP_com' to data source 
  `resul/geom/GP_com.geojson' using driver `GeoJSON'
Writing 150 features with 10 fields and geometry type Polygon.

Importation via une API

Une solution beaucoup plus rapide consiste à importer des données géographiques en faisant appel à une API. Celles-ci renvoient en général des fichiers au format Geojson qu’il sera facile de convertir ensuite au format sf.

A titre d’exemple, nous allons utiliser le site Opendatasoft pour importer des données concernant l’indice de défavorisation sociale à l’échelle des IRIS dans la commune de Fontenay-sous-Bois

Indice de déprivation sociale

Identification du lien de téléchargement

On procède à une sélection des IRIS de la région Ile-de-France puis on copie le lien permettant de récupérer les données au format shapefile

lien de télécargement

Récupération du fichier au format sf

myurl <- "https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/indice-de-defavorisation-sociale-fdep-par-iris/exports/geojson?lang=fr&refine=c_nom_com%3A%22FONTENAY-SOUS-BOIS%22&timezone=Europe%2FBerlin"
map <- st_read(myurl)
Reading layer `OGRGeoJSON' from data source 
  `https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/indice-de-defavorisation-sociale-fdep-par-iris/exports/geojson?lang=fr&refine=c_nom_com%3A%22FONTENAY-SOUS-BOIS%22&timezone=Europe%2FBerlin' 
  using driver `GeoJSON'
Simple feature collection with 21 features and 14 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2.447158 ymin: 48.8391 xmax: 2.500071 ymax: 48.86146
Geodetic CRS:  WGS 84

Visualisation du fonds de carte

plot(map$geometry)

Cartographie rapide d’indicateurs

plot(map["t1_rev_med"], main="Revenu médian")

plot(map["t1_txbac09"], main="Diplômes du supérieur")

plot(map["t1_txchom0"], main="Taux de chômage")

plot(map["t1_txouvr0"], main="Part des ouvriers")

Exercice

  1. Construire un indicateur synthétique combinant les quatre variables
  2. Cartographier cet indicateur
  3. Construire une fonction applicable à une commune quelconque.